library(tibble)
library(forecast)
library(lmtest)

# Panel Values ##
df <- tibble::tibble(
  month_year = as.Date(c("2022-01-01", "2022-02-01", "2022-03-01", "2022-04-01", "2022-05-01", 
                         "2022-06-01", "2022-07-01", "2022-08-01", "2022-09-01", "2022-10-01",
                         "2022-11-01", "2022-12-01", "2023-01-01", "2023-02-01", "2023-03-01", "2023-04-01")),
  Y = c(83.648151, 84.365168, 85.009849, 85.124914, 84.941094, 
        84.902026, 85.067693, 85.459784, 87.319784, 87.412304, 
        86.365590, 85.120584, 84.282976, 84.579300, 84.761890, 83.962681),
  Time = 0:15
)

# Function to run ITS for different intervention start times
run_its <- function(intervention_start) {
  df$Intervention <- ifelse(df$Time >= intervention_start, 1, 0)
  df$TimeAfterIntervention <- ifelse(df$Time >= intervention_start, df$Time - intervention_start, 0)
  
  # Fit ARIMA model with ITS
  model <- lm(Y ~ Time + Intervention + TimeAfterIntervention, data = df)
  
  # Extract coefficients and p-values safely
  model_test <- coeftest(model)
  coeff_names <- rownames(model_test)
  
  # Find correct indices dynamically
  intervention_index <- grep("Intervention", coeff_names)
  time_after_index <- grep("TimeAfterIntervention", coeff_names)
  
  # Safely extract coefficients and p-values
  intervention_coeff <- ifelse(length(intervention_index) > 0, model_test[intervention_index, 1], NA)
  intervention_pval <- ifelse(length(intervention_index) > 0, model_test[intervention_index, 4], NA)
  
  time_after_coeff <- ifelse(length(time_after_index) > 0, model_test[time_after_index, 1], NA)
  time_after_pval <- ifelse(length(time_after_index) > 0, model_test[time_after_index, 4], NA)
  
  return(list(
    Intervention_Coeff = intervention_coeff,
    Intervention_PValue = intervention_pval,
    TimeAfterIntervention_Coeff = time_after_coeff,
    TimeAfterIntervention_PValue = time_after_pval
  ))
}

# Run ITS for intervention shifts (-3 to +3 months)
shifts <- -3:3
results <- lapply(shifts, function(shift) run_its(10 + shift))

# Convert results to a data frame
its_results_df <- tibble::tibble(
  Shift_Months = shifts,
  Intervention_Coeff = sapply(results, function(x) x$Intervention_Coeff),
  Intervention_PValue = sapply(results, function(x) x$Intervention_PValue),
  TimeAfterIntervention_Coeff = sapply(results, function(x) x$TimeAfterIntervention_Coeff),
  TimeAfterIntervention_PValue = sapply(results, function(x) x$TimeAfterIntervention_PValue)
)

# Display results
print(its_results_df)



#########Decahose Values ######

df2 <- tibble::tibble(
  month_year = as.Date(c("2022-01-01", "2022-02-01", "2022-03-01", "2022-04-01", "2022-05-01", 
                         "2022-06-01", "2022-07-01", "2022-08-01", "2022-09-01", "2022-10-01",
                         "2022-11-01", "2022-12-01", "2023-01-01", "2023-02-01", "2023-03-01", "2023-04-01")),
  Y = c(84.6309, 84.6327, 85.6863, 85.6369, 84.9082, 
        85.0195, 85.0075, 85.4599, 85.6809, 85.6093, 
        84.8489, 83.8767, 82.5608, 83.5648, 84.2673, 82.5457),
  Time = 0:15,
  Intervention = c(rep(0, 10), rep(1, 6)), 
  TimeAfterIntervention = c(rep(0, 10), 0:5)
)
its_model2 <- lm(Y ~ Time + Intervention + TimeAfterIntervention, data = df2)
# Function to run ITS for different intervention start times
# Function to run ITS for different intervention start times
run_its <- function(intervention_start) {
  df2$Intervention <- ifelse(df2$Time >= intervention_start, 1, 0)
  df2$TimeAfterIntervention <- ifelse(df2$Time >= intervention_start, df2$Time - intervention_start, 0)
  
  # Fit ARIMA model with ITS
  model2 <- lm(Y ~ Time + Intervention + TimeAfterIntervention, data = df2)
  
  # Extract coefficients and p-values safely
  model_test <- coeftest(model2)
  coeff_names <- rownames(model_test)
  
  # Find correct indices dynamically
  intervention_index <- grep("Intervention", coeff_names)
  time_after_index <- grep("TimeAfterIntervention", coeff_names)
  
  # Safely extract coefficients and p-values
  intervention_coeff <- ifelse(length(intervention_index) > 0, model_test[intervention_index, 1], NA)
  intervention_pval <- ifelse(length(intervention_index) > 0, model_test[intervention_index, 4], NA)
  
  time_after_coeff <- ifelse(length(time_after_index) > 0, model_test[time_after_index, 1], NA)
  time_after_pval <- ifelse(length(time_after_index) > 0, model_test[time_after_index, 4], NA)
  
  return(list(
    Intervention_Coeff = intervention_coeff,
    Intervention_PValue = intervention_pval,
    TimeAfterIntervention_Coeff = time_after_coeff,
    TimeAfterIntervention_PValue = time_after_pval
  ))
}

# Run ITS for intervention shifts (-3 to +3 months)
shifts <- -3:3
results <- lapply(shifts, function(shift) run_its(10 + shift))

# Convert results to a data frame
its_results_df2 <- tibble::tibble(
  Shift_Months = shifts,
  Intervention_Coeff = sapply(results, function(x) x$Intervention_Coeff),
  Intervention_PValue = sapply(results, function(x) x$Intervention_PValue),
  TimeAfterIntervention_Coeff = sapply(results, function(x) x$TimeAfterIntervention_Coeff),
  TimeAfterIntervention_PValue = sapply(results, function(x) x$TimeAfterIntervention_PValue)
)

# Display results
print(its_results_df2)


